Сделайте копию датасета, в которой удалите колонки с количеством пропусков больше 100, а затем удалите все строки с пропусками.
ds <- readRDS('very_low_birthweight.RDS')
glimpse(ds)
## Rows: 671
## Columns: 26
## $ birth <dbl> 81.511, 81.514, 81.552, 81.558, 81.593, 81.602, 81.610, NA, N…
## $ exit <dbl> 81.604, 81.539, 81.552, 81.667, 81.599, 81.771, 81.697, NA, N…
## $ hospstay <int> 34, 9, -2, 40, 2, 62, 32, NA, NA, 28, 38, NA, 62, 69, 1, 93, …
## $ lowph <dbl> NA, 7.250000, 7.059998, 7.250000, 6.969997, 7.189999, 7.32000…
## $ pltct <int> 100, 244, 114, 182, 54, NA, 282, NA, NA, 153, 229, NA, 182, 3…
## $ race <fct> white, white, black, black, black, white, black, NA, NA, blac…
## $ bwt <int> 1250, 1370, 620, 1480, 925, 940, 1255, 600, 700, 1350, 1310, …
## $ gest <dbl> 35.0, 32.0, 23.0, 32.0, 28.0, 28.0, 29.5, 26.0, 24.0, 34.0, 3…
## $ inout <fct> born at Duke, born at Duke, born at Duke, born at Duke, born …
## $ twn <int> 0, 0, 0, 0, 0, 0, 0, NA, NA, 0, 0, NA, 0, 0, 0, 0, 0, 1, 1, 0…
## $ lol <int> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
## $ magsulf <int> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
## $ meth <int> 0, 1, 0, 1, 0, 1, 1, NA, NA, 1, 0, NA, 0, 0, 1, 1, 1, 1, 1, 1…
## $ toc <int> 0, 0, 1, 0, 0, 0, 0, NA, NA, 0, 0, NA, 0, 0, 1, 1, 0, 0, 0, 1…
## $ delivery <fct> abdominal, abdominal, vaginal, vaginal, abdominal, abdominal,…
## $ apg1 <int> 8, 7, 1, 8, 5, 8, 9, NA, NA, 4, 6, NA, 6, 6, 2, 4, 8, 7, 1, 8…
## $ vent <int> 0, 0, 1, 0, 1, 1, 0, NA, NA, 0, 1, NA, 0, 0, 1, 1, 0, 0, 1, 1…
## $ pneumo <int> 0, 0, 0, 0, 1, 0, 0, NA, NA, 0, 0, NA, 1, 0, 1, 0, 0, 0, 1, 1…
## $ pda <int> 0, 0, 0, 0, 0, 0, 0, NA, NA, 0, 0, NA, 0, 0, 0, 0, 0, 0, 0, 0…
## $ cld <int> 0, 0, NA, 0, 0, 0, 0, NA, NA, 0, 0, NA, 1, 0, 0, 1, 0, 0, 0, …
## $ pvh <fct> NA, NA, NA, NA, definite, absent, NA, NA, absent, NA, NA, NA,…
## $ ivh <fct> NA, NA, NA, NA, definite, absent, NA, NA, absent, NA, NA, NA,…
## $ ipe <fct> NA, NA, NA, NA, NA, absent, NA, NA, absent, NA, NA, NA, absen…
## $ year <dbl> 81.51196, 81.51471, 81.55304, 81.55847, 81.59406, 81.60229, 8…
## $ sex <fct> female, female, female, male, female, female, female, NA, NA,…
## $ dead <int> 0, 0, 1, 0, 1, 0, 0, 1, 1, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0…
summary(ds)
## birth exit hospstay lowph
## Min. :81.51 Min. :68.53 Min. :-6574.00 Min. :6.530
## 1st Qu.:83.52 1st Qu.:83.58 1st Qu.: 16.00 1st Qu.:7.130
## Median :84.90 Median :84.96 Median : 37.00 Median :7.210
## Mean :84.75 Mean :84.84 Mean : 40.36 Mean :7.202
## 3rd Qu.:86.07 3rd Qu.:86.17 3rd Qu.: 62.00 3rd Qu.:7.310
## Max. :87.48 Max. :96.87 Max. : 3668.00 Max. :7.550
## NA's :21 NA's :31 NA's :31 NA's :62
## pltct race bwt gest
## Min. : 16.0 black :369 Min. : 400 Min. :22.00
## 1st Qu.:143.0 native American: 16 1st Qu.: 900 1st Qu.:27.00
## Median :202.0 oriental : 4 Median :1120 Median :29.00
## Mean :201.6 white :257 Mean :1094 Mean :28.87
## 3rd Qu.:252.0 NA's : 25 3rd Qu.:1310 3rd Qu.:31.00
## Max. :571.0 Max. :1580 Max. :40.00
## NA's :70 NA's :2 NA's :4
## inout twn lol magsulf
## born at Duke:547 Min. :0.0000 Min. : 0.000 Min. :0.0000
## transported :121 1st Qu.:0.0000 1st Qu.: 0.000 1st Qu.:0.0000
## NA's : 3 Median :0.0000 Median : 3.500 Median :0.0000
## Mean :0.2074 Mean : 8.438 Mean :0.1344
## 3rd Qu.:0.0000 3rd Qu.: 9.000 3rd Qu.:0.0000
## Max. :1.0000 Max. :192.000 Max. :1.0000
## NA's :20 NA's :381 NA's :247
## meth toc delivery apg1
## Min. :0.0000 Min. :0.0000 abdominal:314 Min. :0.000
## 1st Qu.:0.0000 1st Qu.:0.0000 vaginal :335 1st Qu.:2.000
## Median :0.0000 Median :0.0000 NA's : 22 Median :5.000
## Mean :0.4372 Mean :0.2248 Mean :4.903
## 3rd Qu.:1.0000 3rd Qu.:0.0000 3rd Qu.:7.000
## Max. :1.0000 Max. :1.0000 Max. :9.000
## NA's :106 NA's :106 NA's :34
## vent pneumo pda cld
## Min. :0.0000 Min. :0.0000 Min. :0.0000 Min. :0.0000
## 1st Qu.:0.0000 1st Qu.:0.0000 1st Qu.:0.0000 1st Qu.:0.0000
## Median :1.0000 Median :0.0000 Median :0.0000 Median :0.0000
## Mean :0.5803 Mean :0.1969 Mean :0.2087 Mean :0.2694
## 3rd Qu.:1.0000 3rd Qu.:0.0000 3rd Qu.:0.0000 3rd Qu.:1.0000
## Max. :1.0000 Max. :1.0000 Max. :1.0000 Max. :1.0000
## NA's :30 NA's :26 NA's :29 NA's :66
## pvh ivh ipe year sex
## absent :360 absent :442 absent :472 Min. :81.51 female:320
## definite:125 definite: 75 definite: 38 1st Qu.:83.52 male :330
## possible: 41 possible: 10 possible: 17 Median :84.91 NA's : 21
## NA's :145 NA's :144 NA's :144 Mean :84.76
## 3rd Qu.:86.07
## Max. :87.48
## NA's :21
## dead
## Min. :0.0000
## 1st Qu.:0.0000
## Median :0.0000
## Mean :0.2146
## 3rd Qu.:0.0000
## Max. :1.0000
##
factors <- c("race","inout","delivery", "sex", "twn", "vent", "pneumo", "pda", "cld", "dead")
ds_filtered0 <- ds %>%
#удаление колонок где больше 100 пропусков
dplyr::select(where(~sum(is.na(.x)) <= 100)) %>%
#перевод категориальных переменных в факторы
mutate(across(all_of(factors), as.factor)) %>%
#перекодировка факторов
mutate(vent = recode(vent, '0' = "No ventilation", '1' = "Ventilation"),
pneumo = recode(pneumo, '0' = 'No pneumothorax', '1' = 'Pneumothorax'),
pda = recode(pda, '0' = 'No patent ductus arteriosus', '1' = 'Patent ductus arteriosus'),
cld = recode(cld, '0' = 'No oxygen', '1' = 'On oxygen'),
dead = recode(dead, '0' = 'Alive', '1' = 'Dead'),
ID =as.factor(row_number())) %>%
#удаление выбросов
filter(if_any(where(is.numeric),
~ .x >= quantile(.x, 0.25, na.rm = TRUE) - 1.5 * IQR(.x, na.rm = TRUE) &
.x <= quantile(.x, 0.75, na.rm = TRUE) + 1.5 * IQR(.x, na.rm = TRUE)))
#удаление стоблцов с NA
ds_filtered <- ds_filtered0 %>% drop_na
summary(ds_filtered)
## birth exit hospstay lowph
## Min. :81.51 Min. :81.05 Min. :-295.00 Min. :6.530
## 1st Qu.:83.43 1st Qu.:83.56 1st Qu.: 21.00 1st Qu.:7.135
## Median :84.77 Median :84.87 Median : 40.00 Median :7.220
## Mean :84.63 Mean :84.76 Mean : 47.04 Mean :7.215
## 3rd Qu.:85.83 3rd Qu.:85.99 3rd Qu.: 64.00 3rd Qu.:7.320
## Max. :87.48 Max. :87.72 Max. : 797.00 Max. :7.550
##
## pltct race bwt gest
## Min. : 16.0 black :303 Min. : 400 Min. :23.00
## 1st Qu.:148.0 native American: 13 1st Qu.: 960 1st Qu.:28.00
## Median :204.0 oriental : 4 Median :1160 Median :29.00
## Mean :204.5 white :211 Mean :1136 Mean :29.25
## 3rd Qu.:256.0 3rd Qu.:1330 3rd Qu.:31.00
## Max. :571.0 Max. :1500 Max. :36.00
##
## inout twn delivery apg1
## born at Duke:448 0:422 abdominal:262 Min. :0.000
## transported : 83 1:109 vaginal :269 1st Qu.:2.000
## Median :5.000
## Mean :5.024
## 3rd Qu.:7.000
## Max. :9.000
##
## vent pneumo pda
## No ventilation:243 No pneumothorax:438 No patent ductus arteriosus:425
## Ventilation :288 Pneumothorax : 93 Patent ductus arteriosus :106
##
##
##
##
##
## cld year sex dead ID
## No oxygen:393 Min. :81.51 female:264 Alive:467 2 : 1
## On oxygen:138 1st Qu.:83.43 male :267 Dead : 64 4 : 1
## Median :84.77 5 : 1
## Mean :84.63 7 : 1
## 3rd Qu.:85.83 10 : 1
## Max. :87.48 11 : 1
## (Other):525
ds_filtered %>%
pivot_longer(cols = where(is.numeric), names_to = "variable", values_to = "value") %>%
ggplot(aes(x = value, fill = inout)) +
geom_density(alpha = 0.5) +
facet_wrap(~ variable, scales = "free") +
labs(title = "Density Plots of Numeric Variables", x = "Value", y = "Density") +
theme_minimal()
Проведите тест на сравнение значений колонки ‘lowph’ между группами в переменной inout. Вид статистического теста определите самостоятельно. Визуализируйте результат через библиотеку ‘rstatix’. Как бы вы интерпретировали результат, если бы знали, что более низкое значение lowph ассоциировано с более низкой выживаемостью?
t_test_result <- ds_filtered %>%
t_test(lowph ~ inout) %>%
add_significance() %>%
mutate(y.position = max(ds_filtered$lowph, na.rm = TRUE) * 1.05)
t_test_result
## # A tibble: 1 × 10
## .y. group1 group2 n1 n2 statistic df p p.signif y.position
## <chr> <chr> <chr> <int> <int> <dbl> <dbl> <dbl> <chr> <dbl>
## 1 lowph born at … trans… 448 83 5.57 111. 1.77e-7 **** 7.93
ds_filtered %>%
ggplot(aes(x = inout, y = lowph)) +
geom_boxplot() +
stat_pvalue_manual(t_test_result, label = "p") +
labs(x = "Inout Group",
y = "Low pH") +
theme_minimal()
Я выбрала т-тест поскольку довольно много наблюдений, его результаты показывают, что уровень pH значимо отличается у детей родившихся в Дьюке и привезенных из других госпиталей. Это может говорить о том, что транспортировка детей ухудшает их прогноз, возможно им не успевают оказать необходимую помощь.
Сделайте новый датафрейм, в котором оставьте только континуальные или ранговые данные, кроме ‘birth’, ‘year’ и ‘exit’. Сделайте корреляционный анализ этих данных. Постройте два любых типа графиков для визуализации корреляций.
continuous_ranked_vars <- ds_filtered %>%
dplyr::select(where(~ is.numeric(.x) || is.ordered(.x))) %>%
dplyr::select(-birth, -year, -exit)
glimpse(continuous_ranked_vars)
## Rows: 531
## Columns: 6
## $ hospstay <int> 9, 40, 2, 32, 28, 38, 62, 69, 1, 93, 44, 66, 65, 44, 70, 85, …
## $ lowph <dbl> 7.250000, 7.250000, 6.969997, 7.320000, 7.160000, 7.039997, 7…
## $ pltct <int> 244, 182, 54, 282, 153, 229, 182, 361, 378, 255, 186, 260, 18…
## $ bwt <int> 1370, 1480, 925, 1255, 1350, 1310, 1110, 1180, 970, 770, 1490…
## $ gest <dbl> 32.0, 32.0, 28.0, 29.5, 34.0, 32.0, 28.0, 28.0, 28.0, 26.0, 3…
## $ apg1 <int> 7, 8, 5, 9, 4, 6, 6, 6, 2, 4, 8, 1, 8, 5, 9, 9, 9, 6, 2, 1, 6…
cor_matrix <- cor(continuous_ranked_vars, use = "pairwise.complete.obs", method = "spearman")
cor_matrix
## hospstay lowph pltct bwt gest apg1
## hospstay 1.0000000 -0.2337033 -0.15825443 -0.4447105 -0.38014775 -0.1255054
## lowph -0.2337033 1.0000000 0.29144157 0.3234317 0.36957021 0.2730238
## pltct -0.1582544 0.2914416 1.00000000 0.2465972 0.08007581 0.2982137
## bwt -0.4447105 0.3234317 0.24659722 1.0000000 0.71038991 0.3182650
## gest -0.3801477 0.3695702 0.08007581 0.7103899 1.00000000 0.2674875
## apg1 -0.1255054 0.2730238 0.29821365 0.3182650 0.26748752 1.0000000
network_plot(cor_matrix)
chart.Correlation(continuous_ranked_vars, histogram = TRUE, pch = 19)
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
Постройте иерархическую кластеризацию на этом датафрейме.
distance_matrix <- dist(continuous_ranked_vars, method = "euclidean")
hclust_result <- hclust(distance_matrix, method = "ward.D2")
label_colors <- ifelse(ds_filtered$dead == "Alive", "blue", "red")
fviz_dend(hclust_result,
show_labels = TRUE,
label_cols = label_colors,
main = "Hierarchical Clustering Dendrogram",
sub = "", xlab = "", ylab = '')
## Warning: The `<scale>` argument of `guides()` cannot be `FALSE`. Use "none" instead as
## of ggplot2 3.3.4.
## ℹ The deprecated feature was likely used in the factoextra package.
## Please report the issue at <https://github.com/kassambara/factoextra/issues>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
Сделайте одновременный график heatmap и иерархической кластеризации. Интерпретируйте результат.
pheatmap(cor_matrix,
clustering_distance_rows = "euclidean",
clustering_distance_cols = "euclidean",
clustering_method = "ward.D2",
display_numbers = TRUE,
number_format = "%.2f",
main = "Heatmap with Hierarchical Clustering")
По данному графику видно, что числовые переменные можно разделить на 2 группы, которые между собой связаны внутри групп, не считая времени проведенного в госпитале, которая обратно коррелирует со всеми остальными переменными. В первой группе вес при рождении и продолжительность беременности, во второй - pH, значение по шкале Апгар и количество тромбоцитов.
Проведите PCA анализ на этих данных. Проинтерпретируйте результат. Нужно ли применять шкалирование для этих данных перед проведением PCA?
Постройте biplot график для PCA. Раскрасьте его по значению колонки ‘dead’.
Переведите последний график в ‘plotly’. При наведении на точку нужно, чтобы отображалось id пациента.
data_scaled <- scale(continuous_ranked_vars)
pca_result <- PCA(data_scaled, graph = FALSE)
В данном случае данные нужно обязательно шкалировать, так как переменные в разных единицах.
pca_ind <- as.data.frame(pca_result$ind$coord) %>%
mutate(dead = ds_filtered$dead, id = as.character(ds_filtered$ID))
pca_var <- as.data.frame(pca_result$var$coord)
pca_var$Variable <- rownames(pca_var)
p <- ggplot() +
geom_point(data = pca_ind, aes(x = Dim.1, y = Dim.2, color = dead, text = id)) +
geom_segment(data = pca_var, aes(x = 0, y = 0, xend = Dim.1, yend = Dim.2)) +
geom_text(data = pca_var, aes(x = Dim.1, y = Dim.2, label = Variable), hjust = 1.2, vjust = 1.2) +
scale_color_manual(values = c("blue", "red")) +
labs(title = "Interactive PCA Biplot", x = "PC1", y = "PC2") +
theme_minimal()
## Warning in geom_point(data = pca_ind, aes(x = Dim.1, y = Dim.2, color = dead, :
## Ignoring unknown aesthetics: text
ggplotly(p, tooltip = "text")
сложно интерпретировать этот график, я пыталась по разному менять размеры стрелок и точек, но для этих данных проще построить отдельный график со стрелками чтобы посмотреть на вклад переменных в главные компоненты
fviz_pca_var(pca_result, col.var = "contrib", gradient.cols = c("blue", "yellow", "red"),
repel = TRUE, title = "Contribution of Variables to PCA")
На этом графике еще более четко видно разделение переменных по 3 группам описанным выше. Видно, что есть 2 переменные - вес и срок беременности, которые дают самый весомый вклад в первую компоненту и она объясняет довольно много дисперсии (почти 40%). Вероятно, это 2 главных фактора, которые влияют на состояние недоношенных младенцев. Другие 3 переменные, которые харакетризиют жизненные показатели младенцев тоже вносят большой вклад в первую компоненту, возможно они свзяны с осложнениями при беременности и родах, таких как преэклапсия или инфекции. Вторая компонента основана главным образом на переменной время проведенное в госпитале, вторая компонента объясняет мало дисперсии данных по сравнению с первой. Возможно потому что время пребывания не всегда связано с состоянием пациента, или эта связь сложная, к примеру мало прибывают те, кто умирает или наоборот в хорошем состоянии, а те, кому нужна помощь будут находиться дольше. Похоже на то, что первая компонента разделяет младенцев на умерших и выживших.
Некорректно использовать колонку dead, так как она не дает информации о том когда именно была смерть.
explained_variance <- pca_result$eig[, 2]
explained_variance_df <- data.frame(
Component = paste0("PC", 1:length(explained_variance)),
Variance = explained_variance
)
ggplot(explained_variance_df, aes(x = Component, y = Variance)) +
geom_bar(stat = "identity", fill = "steelblue") +
geom_text(aes(label = paste0(round(Variance, 1), "%")), vjust = -0.5, size = 4) +
labs(
title = "Explained Variance by Principal Components",
x = "Principal Component",
y = "Explained Variance (%)"
) +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
Еще хороший график чтобы понять какие кроме первых двух компонент могут еще вносить вклад. В нашем случае не супер информативно, тк мало переменных в исходном датасете.
umap_result <- umap(continuous_ranked_vars, n_neighbors = 15, min_dist = 0.1, n_components = 2)
umap_data <- as.data.frame(umap_result$layout)
umap_data$dead <- ds_filtered$dead
ggplot(umap_data, aes(x = V1, y = V2, color = dead)) +
geom_point(size = 3, alpha = 0.7) +
labs(title = "UMAP Projection", x = "UMAP 1", y = "UMAP 2") +
theme_minimal()
С помощью UMAP данные лучше разделились по 2 компоненте, чем по первой, умершие попали так же преимущественно в один кластер как и на PCA.
Чтобы проверить как влияют параметры я сделала 4 варианта параметров и запустила для них UMAP
umap_configs <- list(
list(n_neighbors = 5, min_dist = 0.1),
list(n_neighbors = 15, min_dist = 0.1),
list(n_neighbors = 15, min_dist = 0.5),
list(n_neighbors = 30, min_dist = 0.5)
)
apply_umap <- function(config, data) {
result <- umap(data, n_neighbors = config$n_neighbors, min_dist = config$min_dist, n_components = 2)
as.data.frame(result$layout) %>%
mutate(n_neighbors = config$n_neighbors, min_dist = config$min_dist)
}
umap_results <- bind_rows(lapply(umap_configs, apply_umap, data = continuous_ranked_vars))
umap_results <- umap_results %>%
mutate(dead = rep(ds_filtered$dead, times = length(umap_configs)))
ggplot(umap_results, aes(x = V1, y = V2, color = dead)) +
geom_point(size = 2, alpha = 0.7) +
facet_wrap(~ n_neighbors + min_dist, nrow = 2, ncol = 2, labeller = label_both) +
labs(title = "UMAP Projections with Different Parameters",
x = "UMAP 1", y = "UMAP 2") +
scale_color_manual(values = c("blue", "red")) +
theme_minimal()
Проекции с разными параметрами имеют разное количество кластеров: чем больше оба параметра - тем меньше кластеров и тем лучше видно как точки внутри кластеров расположены относительно друг друга. На самом первом графике максимальное количество групп точек, так как алгоритм учитывает только 5 соседних наблюдений, на самом же последнем графике точки все точки сгруппированы вместе.
Пермутация
set.seed(123)
data_perm_50 <- continuous_ranked_vars %>%
mutate(bwt = ifelse(runif(n()) <= 0.5, sample(bwt), bwt))
data_perm_100 <- continuous_ranked_vars %>%
mutate(bwt = sample(bwt))
Функции для PCA и UMAP
run_pca <- function(data) {
pca_result <- PCA(data, graph = FALSE)
return(pca_result)
}
run_umap <- function(data) {
umap_result <- umap(data, n_neighbors = 15, min_dist = 0.1, n_components = 2)
return(as.data.frame(umap_result$layout))
}
#PCA и UMAP для оригинальных данных
pca_original <- run_pca(continuous_ranked_vars)
umap_original <- run_umap(continuous_ranked_vars)
#PCA и UMAP для данных с пермутацией 50%
pca_perm_50 <- run_pca(data_perm_50)
umap_perm_50 <- run_umap(data_perm_50)
#PCA и UMAP для данных с пермутацией 100%
pca_perm_100 <- run_pca(data_perm_100)
umap_perm_100 <- run_umap(data_perm_100)
explained_variance <- tibble(
Component = paste0("PC", 1:length(pca_original$eig[, 2])),
Original = pca_original$eig[, 2],
Permuted_50 = pca_perm_50$eig[, 2] ,
Permuted_100 = pca_perm_100$eig[, 2]
)
print(explained_variance)
## # A tibble: 6 × 4
## Component Original Permuted_50 Permuted_100
## <chr> <dbl> <dbl> <dbl>
## 1 PC1 39.4 32.7 30.3
## 2 PC2 17.9 17.4 18.0
## 3 PC3 15.0 15.3 16.3
## 4 PC4 12.2 14.1 14.3
## 5 PC5 11.0 12.0 12.0
## 6 PC6 4.47 8.45 9.03
Из таблицы выше видно, что пермутация данных ухудшила PCA - первая компонента стала объяснять меньше вариабельности почти на 25% для 100% пермутированных данных по сравнению с оригинальными.
pca_data <- list(
Original = as.data.frame(pca_original$ind$coord) %>% mutate(Method = "PCA", Data = "Original", dead = ds_filtered$dead),
Permuted_50 = as.data.frame(pca_perm_50$ind$coord) %>% mutate(Method = "PCA", Data = "Permuted 50%", dead = ds_filtered$dead),
Permuted_100 = as.data.frame(pca_perm_100$ind$coord) %>% mutate(Method = "PCA", Data = "Permuted 100%", dead = ds_filtered$dead)
) %>% bind_rows()
umap_data <- list(
Original = umap_original %>%
mutate(Method = "UMAP", Data = "Original", dead = ds_filtered$dead),
Permuted_50 = umap_perm_50 %>%
mutate(Method = "UMAP", Data = "Permuted 50%", dead = ds_filtered$dead),
Permuted_100 = umap_perm_100 %>%
mutate(Method = "UMAP", Data = "Permuted 100%", dead = ds_filtered$dead)
) %>% bind_rows()
combined_data <- bind_rows(
pca_data %>% rename(Dim1 = Dim.1, Dim2 = Dim.2),
umap_data %>% rename(Dim1 = V1, Dim2 = V2)
)
ggplot(combined_data, aes(x = Dim1, y = Dim2, color = dead)) +
geom_point(size = 2, alpha = 0.5) +
facet_grid(Method ~ Data) +
labs(title = "PCA and UMAP Projections with Permuted Data",
x = "Dimension 1", y = "Dimension 2") +
theme_minimal() +
theme(legend.position = "bottom")
Честно говоря на этих данных сложно увидеть эффект пермутации, в теории кластеры должны были значительно поменяться, но на PCA эффект не очень заметен, на UMAP проекции с пермутациями отличаются от изначальной тем, что точки для мертвых не группируются вместе. Чтобы лучше увидеть эффект для PCA нарисовала еще графики со стрелками, на них видно, что bwt вносит меньше вклада чем больше пермутированных значений.
pca_var_original <- fviz_pca_var(pca_original,
col.var = "contrib",
gradient.cols = c("blue", "green", "red"),
repel = TRUE,
title = "PCA Variables: Original Data")
# Permuted 50% Data
pca_var_perm_50 <- fviz_pca_var(pca_perm_50,
col.var = "contrib",
gradient.cols = c("blue", "green", "red"),
repel = TRUE,
title = "PCA Variables: Permuted 50%")
# Permuted 100% Data
pca_var_perm_100 <- fviz_pca_var(pca_perm_100,
col.var = "contrib",
gradient.cols = c("blue", "green", "red"),
repel = TRUE,
title = "PCA Variables: Permuted 100%")
grid.arrange(pca_var_original, pca_var_perm_50, pca_var_perm_100, nrow = 1)
Давайте проведем анализ чувствительности. Проведите анализ, как в шагах 4-6 для оригинального с удалением всех строк с пустыми значениями (т.е. включая колонки с количеством пропущенных значений больше 100), а затем для оригинального датафрейма с импутированием пустых значений средним или медианой. Как отличаются получившиеся результаты? В чем преимущества и недостатки каждого подхода?
#импутировала NA регрессией с учетом других значений
vars_to_exclude <- c("ID", "year", "birth", "exit")
ds_filtered_imp0 <- ds_filtered0 %>%
dplyr::select(-all_of(vars_to_exclude)) %>%
mice() %>%
complete()
##
## iter imp variable
## 1 1 hospstay lowph pltct race bwt gest inout twn delivery apg1 vent pneumo pda cld sex
## 1 2 hospstay lowph pltct race bwt gest inout twn delivery apg1 vent pneumo pda cld sex
## 1 3 hospstay lowph pltct race bwt gest inout twn delivery apg1 vent pneumo pda cld sex
## 1 4 hospstay lowph pltct race bwt gest inout twn delivery apg1 vent pneumo pda cld sex
## 1 5 hospstay lowph pltct race bwt gest inout twn delivery apg1 vent pneumo pda cld sex
## 2 1 hospstay lowph pltct race bwt gest inout twn delivery apg1 vent pneumo pda cld sex
## 2 2 hospstay lowph pltct race bwt gest inout twn delivery apg1 vent pneumo pda cld sex
## 2 3 hospstay lowph pltct race bwt gest inout twn delivery apg1 vent pneumo pda cld sex
## 2 4 hospstay lowph pltct race bwt gest inout twn delivery apg1 vent pneumo pda cld sex
## 2 5 hospstay lowph pltct race bwt gest inout twn delivery apg1 vent pneumo pda cld sex
## 3 1 hospstay lowph pltct race bwt gest inout twn delivery apg1 vent pneumo pda cld sex
## 3 2 hospstay lowph pltct race bwt gest inout twn delivery apg1 vent pneumo pda cld sex
## 3 3 hospstay lowph pltct race bwt gest inout twn delivery apg1 vent pneumo pda cld sex
## 3 4 hospstay lowph pltct race bwt gest inout twn delivery apg1 vent pneumo pda cld sex
## 3 5 hospstay lowph pltct race bwt gest inout twn delivery apg1 vent pneumo pda cld sex
## 4 1 hospstay lowph pltct race bwt gest inout twn delivery apg1 vent pneumo pda cld sex
## 4 2 hospstay lowph pltct race bwt gest inout twn delivery apg1 vent pneumo pda cld sex
## 4 3 hospstay lowph pltct race bwt gest inout twn delivery apg1 vent pneumo pda cld sex
## 4 4 hospstay lowph pltct race bwt gest inout twn delivery apg1 vent pneumo pda cld sex
## 4 5 hospstay lowph pltct race bwt gest inout twn delivery apg1 vent pneumo pda cld sex
## 5 1 hospstay lowph pltct race bwt gest inout twn delivery apg1 vent pneumo pda cld sex
## 5 2 hospstay lowph pltct race bwt gest inout twn delivery apg1 vent pneumo pda cld sex
## 5 3 hospstay lowph pltct race bwt gest inout twn delivery apg1 vent pneumo pda cld sex
## 5 4 hospstay lowph pltct race bwt gest inout twn delivery apg1 vent pneumo pda cld sex
## 5 5 hospstay lowph pltct race bwt gest inout twn delivery apg1 vent pneumo pda cld sex
ds_filtered_imp <- bind_cols(ds_filtered_imp0, ds_filtered0 %>% dplyr::select(all_of(vars_to_exclude)))
continuous_ranked_vars_imp <- ds_filtered_imp %>%
dplyr::select(where(~ is.numeric(.x) || is.ordered(.x))) %>%
dplyr::select(-birth, -year, -exit)
cor_matrix <- cor(continuous_ranked_vars_imp, use = "pairwise.complete.obs", method = "spearman")
network_plot(cor_matrix)
distance_matrix <- dist(continuous_ranked_vars_imp, method = "euclidean")
hclust_result <- hclust(distance_matrix, method = "ward.D2")
label_colors <- ifelse(ds_filtered0$dead == "Alive", "blue", "red")
fviz_dend(hclust_result,
show_labels = TRUE,
label_cols = label_colors,
main = "Hierarchical Clustering Dendrogram",
sub = "", xlab = "", ylab = '')
pheatmap(cor_matrix,
clustering_distance_rows = "euclidean",
clustering_distance_cols = "euclidean",
clustering_method = "ward.D2",
display_numbers = TRUE,
number_format = "%.2f",
main = "Heatmap with Hierarchical Clustering")
chart.Correlation(continuous_ranked_vars_imp, histogram = TRUE, pch = 19)
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
chart.Correlation(continuous_ranked_vars, histogram = TRUE, pch = 19)
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
Лучше всего видно как повлияла импутация на двух последних графиках: значения корреляции немного поменялись, для некоторых пар значимная корреляция есть только на одном из графиков. Наибольшие изменения для переменных gest и hospstay, у них даже форма гистограмм поменялась. Но в данном анализе сложно понять что больше повлияло - импутация или то что больше 100 строк данных добавилось.
data_scaled <- scale(continuous_ranked_vars_imp)
pca_result <- PCA(data_scaled, graph = FALSE)
fviz_pca_var(pca_result, col.var = "contrib", gradient.cols = c("blue", "yellow", "red"),
repel = TRUE, title = "Contribution of Variables to PCA")
pca_ind <- as.data.frame(pca_result$ind$coord) %>%
mutate(dead = ds_filtered0$dead, id = as.character(ds_filtered0$ID))
pca_var <- as.data.frame(pca_result$var$coord)
pca_var$Variable <- rownames(pca_var)
p <- ggplot() +
geom_point(data = pca_ind, aes(x = Dim.1, y = Dim.2, color = dead, text = id)) +
geom_segment(data = pca_var, aes(x = 0, y = 0, xend = Dim.1, yend = Dim.2)) +
geom_text(data = pca_var, aes(x = Dim.1, y = Dim.2, label = Variable), hjust = 1.2, vjust = 1.2) +
scale_color_manual(values = c("blue", "red")) +
labs(title = "Interactive PCA Biplot", x = "PC1", y = "PC2") +
theme_minimal()
## Warning in geom_point(data = pca_ind, aes(x = Dim.1, y = Dim.2, color = dead, :
## Ignoring unknown aesthetics: text
ggplotly(p, tooltip = "text")
umap_result <- umap(continuous_ranked_vars_imp, n_neighbors = 15, min_dist = 0.1, n_components = 2)
umap_data <- as.data.frame(umap_result$layout)
umap_data$dead <- ds_filtered0$dead
ggplot(umap_data, aes(x = V1, y = V2, color = dead)) +
geom_point(size = 3, alpha = 0.7) +
labs(title = "UMAP Projection", x = "UMAP 1", y = "UMAP 2") +
theme_minimal()
Видно, что после импутации первая компонента PCA стала объяснять больше изменчивости и на графике проекции умершие и выжившие разделяются на группы гораздо лучше, что говорит о том, что по этим данным лучше получится определить связь смертности с указаными переменными. Еще интересно, что появились достаточно четкие выбросы на PCA, но их нет на UMAP. На UMAP по ощущению импутация вообще не повлияла.
Интересно, что в лекциях не было ничего сказано, что t-SNE, который похож на UMAP, но в биоинформатике гораздо более популярен. Попробовала им тоже сравнить данные до и после импутации и получилось довольно сильное отличие.
set.seed(42)
tsne_result <- Rtsne(as.matrix(continuous_ranked_vars), dims = 2, perplexity = 30, verbose = TRUE, max_iter = 500)
## Performing PCA
## Read the 531 x 6 data matrix successfully!
## OpenMP is working. 1 threads.
## Using no_dims = 2, perplexity = 30.000000, and theta = 0.500000
## Computing input similarities...
## Building tree...
## Done in 0.06 seconds (sparsity = 0.209277)!
## Learning embedding...
## Iteration 50: error is 56.340360 (50 iterations in 0.03 seconds)
## Iteration 100: error is 53.994127 (50 iterations in 0.03 seconds)
## Iteration 150: error is 53.890921 (50 iterations in 0.03 seconds)
## Iteration 200: error is 53.880156 (50 iterations in 0.03 seconds)
## Iteration 250: error is 53.876007 (50 iterations in 0.03 seconds)
## Iteration 300: error is 0.592579 (50 iterations in 0.03 seconds)
## Iteration 350: error is 0.507048 (50 iterations in 0.03 seconds)
## Iteration 400: error is 0.491431 (50 iterations in 0.03 seconds)
## Iteration 450: error is 0.482132 (50 iterations in 0.03 seconds)
## Iteration 500: error is 0.471446 (50 iterations in 0.03 seconds)
## Fitting performed in 0.29 seconds.
tsne_data <- data.frame(tsne_result$Y, dead = ds_filtered$dead)
colnames(tsne_data) <- c("Dim1", "Dim2", "dead")
# t-SNE визуализация
ggplot(tsne_data, aes(x = Dim1, y = Dim2, color = dead)) +
geom_point(size = 2, alpha = 0.7) +
labs(title = "t-SNE Visualization", x = "Dimension 1", y = "Dimension 2") +
theme_minimal()
set.seed(42)
tsne_result <- Rtsne(as.matrix(continuous_ranked_vars_imp), dims = 2, perplexity = 30, verbose = TRUE, max_iter = 500)
## Performing PCA
## Read the 670 x 6 data matrix successfully!
## OpenMP is working. 1 threads.
## Using no_dims = 2, perplexity = 30.000000, and theta = 0.500000
## Computing input similarities...
## Building tree...
## Done in 0.09 seconds (sparsity = 0.166180)!
## Learning embedding...
## Iteration 50: error is 59.616393 (50 iterations in 0.06 seconds)
## Iteration 100: error is 55.926854 (50 iterations in 0.04 seconds)
## Iteration 150: error is 55.684807 (50 iterations in 0.04 seconds)
## Iteration 200: error is 55.658898 (50 iterations in 0.04 seconds)
## Iteration 250: error is 55.665700 (50 iterations in 0.04 seconds)
## Iteration 300: error is 0.714897 (50 iterations in 0.04 seconds)
## Iteration 350: error is 0.582656 (50 iterations in 0.04 seconds)
## Iteration 400: error is 0.546329 (50 iterations in 0.05 seconds)
## Iteration 450: error is 0.534019 (50 iterations in 0.05 seconds)
## Iteration 500: error is 0.523819 (50 iterations in 0.05 seconds)
## Fitting performed in 0.45 seconds.
tsne_data <- data.frame(tsne_result$Y, dead = ds_filtered0$dead)
colnames(tsne_data) <- c("Dim1", "Dim2", "dead")
# t-SNE визуализация
ggplot(tsne_data, aes(x = Dim1, y = Dim2, color = dead)) +
geom_point(size = 2, alpha = 0.7) +
labs(title = "t-SNE Visualization", x = "Dimension 1", y = "Dimension 2") +
theme_minimal()